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A structure prediction method is presented based on the Minima Hopping method. Optimized moves on the 
configurational enthalpy surface are performed to escape local minima using variable cell shape molecular 
dynamics by aligning the initial atomic and cell velocities to low curvature directions of the current mini- 
mum. The method is applied to both silicon crystals and binary Lennard- Jones mixtures and the results are 
compared to previous investigations. It is shown that a high success rate is achieved and a reliable prediction 
of unknown ground state structures is possible. 
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I. INTRODUCTION 

Predicting structures of periodic systems^ is one of the 
most challenging tasks in material sciences. Not only in 
material design under various external conditions, but 
also in biology and pharmacy there is an increasing de- 
mand for efficient and reliable structure prediction meth- 
ods»2i£ A very common way of predicting favorable struc- 
tures is extracting known structures from databases of 
structures previously found in similar materials. The en- 
ergetically most stable structure is identified and gives 
the putative ground state. However, this approach has 
a limited success rate when the ground state is an un- 
known structure, which can only be found by performing 
an extensive search. Similarly, data mining is capable of 
predicting new crystalline structures based on a huge set 
of experimental data and/or ab initio calculations!^ 

Recently, more advanced methods for crystal structure 
prediction have been developed and applied, which allow 
a systematic search for the ground state structure based 
solely on the system's composition and the external con- 
ditions. The most promising of these methods and their 
applications on crystal structure prediction include simu- 
lated annealingj£~— genetic algorithms^— and metady- 
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Some of the methods, their advantages and 
drawbacks are characterized and discussed by Oganov 
et al.;& giving a good overview over the aforementioned 
structure prediction methods. 

The Minima Hopping (MH) method is an algorithm 
which allows an efficient exploration of a high dimen- 
sional potential energy surface of complex systems, while 
progressing toward the global minimum structured It 
has been successfully applied to various isolated systems 
such as Lennard- Jones and silicon clusters r^— doped 
silicon fullerenes^ complex biological molecules^ and 
large gold clusters. 31 The MH method has also been used 
in other applications, such as providing realistic atomic 
force microscopy (AFM) tips for AFM simulations^ - — 
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In this paper we present an approach for structure pre- 
diction of crystals by generalizing the MH method to pe- 
riodic systems with pressure constraint. The method was 
tested on two benchmark systems, silicon crystals and bi- 
nary Lennard- Jones (BLJ) mixtures. Finding the ground 
state of the Si64 cell presents a challenging task due to 
its multi-funnel character. Volume constraints are used 
to demonstrate the capability of predicting porous silicon 
crystal structures. Softening, a procedure to modify ini- 
tial molecular dynamics velocities, is used to increase the 
efficiency. To demonstrate the effect of softening statisti- 
cal data are collected for its application on a BLJ mixture 
in a small cell. Furthermore, a larger, much-studied BLJ 
mixture is investigated and a new putative ground state 
structure is presented. 



II. METHOD AND APPLICATION 

A. Generalizing the Minima Hopping to 
periodic systems 

The initial implementation of the MH method was 
capable of performing a search for low-lying min- 
ima on the SA^-dimensional potential energy surface 
E = E{ri), i = 1, . . . , N, of an isolated molecule or a pe- 
riodic system in a rigid box with N atoms. Starting from 
some current local minimum on the potential energy sur- 
face, an escape step is performed using a short molecular 
dynamics (MD) simulation which is stopped as soon as 
mdmin potential energy minima have been found along 
the MD trajectory. Then, a local geometry relaxation 
is performed. The new minimum is either accepted or 
rejected, depending on its energy and a threshold value 
Ediff, which is updated by a feedback such that half of 
the new configurations are accepted. To avoid exploring 
previously visited regions of the potential energy surface, 
the kinetic energy Ekm is increased when already known 
structures are revisited. Then, the cycle is repeated start- 
ing with a new MD escape. 

However, it is essential for the purpose of structure 
prediction in periodic systems that not only the atomic 



positions are allowed to be optimized, but also the cell 
shape, especially when external constraints are imposed. 
Furthermore, using a variable cell shape reduces the en- 
thalpy barriers for phase transitions. Hence, to gener- 
alize the MH method for periodic systems with variable 
cell shape, the degrees of freedom are augmented by the 
three variable cell vectors a, b and c, describing the edges 
of the simulation cell. They are combined to a 3 x 3 
tensor h = {a, b, c}. The atomic positions can be ex- 
pressed by vectors in lattice coordinates Sj according to 
r.j = h&i- The potential energy surface E — E(s il h) is 
now a function of the cell vectors and the atomic posi- 
tions with respect to the cell vectors. When an external 
pressure P is applied, the potential energy is replaced 
by the conhgurational enthalpy H c — E(si, h) + Pfl(h), 
where $7 = det(h) is the volume of the simulation cell. 
Periodic boundary conditions are applied for the evalua- 
tion of H c . 

Finding the local extrema of the conhgurational en- 
thalpy is equivalent to finding the set of coordinates for 
which the following conditions are satisfied: 
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where sj, 7 € {x,y,z} are the components of the frac- 
tional coordinate of atom i and h a p denote the elements 
of the tensor h. 

Similarly, for the escape step, the MD needs to be 
performed taking into account the additional cell pa- 
rameters. Hence, both the atomic positions in lattice 
coordinates s,(t) and the lattice vectors h(t) are time- 
dependent. A Lagrangian to perform variable cell shape 
MD at constant pressure P was proposed by Parrincllo 
and Rahman in 1980 35 and has been widely applied: 
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where g — h T h is the metric tensor, W is the fictitious 
mass of the simulation cell and rrii is the mass of atom 
i. The first two terms in equation ([3J correspond to the 
kinetic and potential energy of the atoms, respectively. 
The third and the fourth term correspond to the kinetic 
and potential energy of the simulation cell, respectively. 
The interaction between the atoms are described by a ra- 



equation (QJ is related to the components of the nega- 
tive volume-normalized strain derivatives of the energy 
a a /3 — —fi-g^r- (the stress tensor for the static case) by 

f h h T = —fla, where e is the strain tensor (see Ref. 36]). 
Using the Lagrange's equation, the equations of motion 
for the atomic positions in lattice coordinates and the cell 
vectors are then given by 
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The tensor in equation (|SJ) can be expanded to 
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Furthermore we have the following relation: 

r| = {bxc,cxa,axb} 



(7) 



(8) 



Since the forces are velocity-dependent, the time in- 
tegration should be performed by an advanced integra- 
tion scheme to correctly find the numerical solution of 
the equations of motion ((3J and © (for example the 
symplectic Runge-Kutta method or a predictor-corrector 
scheme). However, we are only interested in escaping the 
local minimum and we only perform short MD runs with 
few oscillations in the potential energy. Therefore, the 
long-time conservation of the total energy is not an issue 
and we can simply use finite differences to discretize the 
equations of motion, equivalent to the Verlet algorithm. 

If we define the following generalized forces 
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and use finite difference formulae for Sj, S{ and h, the time 
evolution of the atomic positions and the cell vectors can 
be computed iteratively according to 



Sl (t + At) = Si (t) + AtVi(t)+ 
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h(t + At) = h(t) + AtV(t) + At 2 —F{t) (12) 
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dial pair potential (/>(r^) where r.y = I r.; — iv,- 1 . However, Here we introduced the dummy variables 
the formalism can readily be extended to any many-body 
potential which are expressed by atomic positions. 
The symmetrized stress tensor n can be written as 



Vift) := «W-ff-^ and V(t) := h(t) - 
which can be considered as velocities of the correspond- 
ing variables. 
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B. Softening and optimizing cell parameters 



where we use the cell gradients of the potential energy 
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and define vw = hs;. The second term in 



Softening, a method of biasing the initial velocities for 
the MD simulation of the escape step, is used to increase 



the efficiency of MH^ 7 - For the crystal structure predic- 
tion MH the velocity vector consists not only of atomic 
velocities, but also of the cell velocities. First, a ran- 
dom velocity direction with Gaussian distributed magni- 
tudes is chosen. The initial velocity amplitude is chosen 
that the kinetic energy is small, allowing only low bar- 
riers to be crossed during the MD escape. In chemistry, 
the Bell-Evans-Polanyi 38 principle states that strongly 
exothermic reactions have a low activation energy. Re- 
actant and product are in fact neighboring local minima 
on the potential energy surface and the chemical reac- 
tion is a transition of an energy barrier connecting the 
two minima. Hence, the Bell-Evans-Polanyi principle can 
be generalized to any transitions between local minima 
on the potential energy surface during MD simulations. 
Recently, Roy et al. have shown that, on average, cross- 
ing low energy barriers along MD trajectories will lead 
into the basin of attraction of lower energy local min- 
ima than crossing high energy barriers. 39 Furthermore, 
low energy barriers are generally connected to low fre- 
quency eigenmodes of local minimal These properties 
can be readily extended to the configurational enthalpy. 
Therefore, the probability of finding low enthalpy config- 
urations can be expected to increase when the direction 
of the initial velocity vector of a MD run points toward 
a direction with low curvature. Hence, in a second step, 
the velocity vector from the first step is rotated such 
that it is oriented along soft mode directions of the cur- 
rent minimum. The rotation procedure is performed by 
iteratively minimizing the energy along the escape di- 
rection at a constant distance from the local minimum. 37 
However, over-biasing the velocities is not favorable since 
the random and therefore ergodic character of the escape 
step should be retained in order not to arrive at a deter- 
ministic process. In fact, softening is mainly applied to 
eliminate components of the velocity on hard modes. 

If a quasi-Newton method (like the Broyden-Fletcher- 
Goldfarb-Shanno method) is used to perform local ge- 
ometry optimization within the MH method an alterna- 
tive softening procedure is possible. In a quasi-Newton 
method an approximate Hessian matrix is continuously 
updated during a geometry relaxation. Before perform- 
ing the MD escape step the approximate Hessian ma- 
trix from the previous relaxation step is diagonalized 
and a small number of low frequency eigenvectors are 
extracted. A randomized superposition of these eigen- 
vectors are then used to provide the initial, soft veloc- 
ity vector for the following MD trajectory. If feasible, 
the Hessian matrix can also be computed analytically to 
evaluate the soft mode directions. 

Softening has been successfully used in previous ap- 
plications of MHi 30 ' 31 ' 37 In Table [I] the performance of 
the MH method with and without softening is compared 
for a benchmark system, a BLJ mixture with type A 
and B atoms in a small cell, AgB^ It can be clearly 
seen that the curvature of the configurational enthalpy 
along the velocity direction is reduced by roughly one 
order of magnitude when softening is used. The sig- 



nificant increase in the efficiency when softening is ap- 
plied can also be ascribed to improving the cell veloc- 
ities, since the impact of cell parameters on the struc- 
ture can be larger than the atomic coordinates. For 
cells with a small number of atoms this can be illus- 
trated by a simple example. Consider a simulation cell 
containing one single atom at the origin. The lattice 
vectors are given by a = |(y + z — x), b = |(z + x — y) 
and c = ^(x + y — z) (where the hats denote the unit 
vectors and a is the lattice constant), which define a 
body-centered lattice. Assume the atomic coordinates 
x were fixed. Transforming the cell to a* = |(t/ + z), 
b* = ~(z + x) and c* = | (x + y) will result in a face- 
centered cubic lattice, a totally different structure. How- 
ever, when the cell parameters a, b and c are fixed, 
there is no possibility to transform the atomic coordi- 
nates x to a system where a face-centered cubic lattice 
is obtained. Obviously, the impact of the cell parame- 
ters decreases with increasing number of atoms, which 
is equivalent to the limit where an infinitely large cell 
is used. Similarly, this principle is also used in metha- 
dynamics simulations where the cell parameters are cho- 
sen as the collective variables^ Methadynamics simula- 
tions were recently applied to predict phase transitions 
in silicon based on neural-network representation of the 
density functional theory potential energy surface, where 
the Gibbs free energy surface was explored at fixed pres- 
sure and temperature as a function of the cell parameters 
h£L The effect of softening on the efficiency of the MH 
method is in general larger for small periodic cells with 
a low number of atoms than for isolated molecules of 
similar size. So, for non-periodic systems with the same 
number of atoms as our benchmark cell softening does 
not give a significant performance improvement. 



The fictitious cell mass W is another adjustable pa- 
rameter. Choosing a too large ratio for cell and atomic 
mass will result in a very stiff cell which will not adjust it- 
self smoothly during the simulation. However, when this 
ratio is to small the cell can fluctuate violently resulting 
in a strong deformation of the cell within one step of the 
simulation. We have found that choosing the cell mass 
similar to the atomic masses is a reasonable choice. 



During a MH simulation it can happen that the cell 
shape gets heavily distorted leading to small angles be- 
tween the three lattice vectors, and thus resulting in a 
very flat cell. This behaviour is not desired since more 
periodic cells have to be taken into account when com- 
puting the potential energy and its derivatives. Further- 
more, it makes it difficult to identify equivalent struc- 
tures in different cells. Therefore, whenever necessary, a 
transformation of the cell vectors is performed to obtain 
shorter cell vectors (for details see Ref. [25]). 



TABLE I. The impact of softening on various quantities of 
the MH method is shown. Starting from 100 different ran- 
dom input configurations all runs were continued until the 
ground state structure was found, resulting in a success rate 
of 100%. The first column shows the different degrees of free- 
dom (DOF) that are taken into account during softening. The 
second and third column show the median value of the curva- 
ture of the enthalpy k along the initial MD velocity direction 
before (Kb) and after (« a ) softening, respectively. The fourth 
column contains the median value of the number of visited 
minima n m i n before reaching the global minimum. A BLJ 
mixture with A8B4 atoms was used described by the modi- 
fied Lennard- Jones potential as discussed in section Hi D I The 
parameters were set to oaa = 1.50, oab = 2.25, ubb = 3.00 
and tAA — 1.00, eab = 1.25, ess = 1.00, and a cutoff radius 
°f r< aJ3 = 2.5<T a £ was used. 

§The number of softening iterations was doubled. 
^Initial atomic velocities were set to zero before softening. 
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C. Application on silicon with optimized 
performance by parallelization and lattice vector 
prediction 

As an application of the MH method, a silicon su- 
percell with 64 atoms was investigated at zero pressure. 
Since the number of local minima increases exponentially 
with respect to the number of atoms in a system, find- 
ing the global minimum of a cell with many atoms is 
a challenging task. Furthermore, the silicon supercell 
with 64 atoms is a multi-funnel system with crystalline 
and amorphous structures and therefore presents an ad- 
ditional challenge for global optimization. For a silicon 
supercell with 64 atoms in a rigid bo:»2& the difficulty of 
finding the ground state has already been demonstrated. 
Several million minima had to be visited before finding 
the ground state, the well-known cubic diamond struc- 
ture. To demonstrate the advantage of the variable cell 
shape MH, we revisit the same problem set. Using the 
environment-dependent interatomic potential (EDIP) for 
silicon^— statistical data were collected for a set of 100 
serial MH runs at zero pressure. Each run started with a 
highly random configuration in a distorted cell and was 
stopped as soon as a structure with ground state energy 
was found or 8000 distinct minima were accepted, a small 
number for such a large system. Since some of the local 
minima are visited several times and only half of the new 
structures are accepted the search length corresponds to 
visiting at most some 50 000 minima. The success rate 
was found to be 45%. Each successful run visited an av- 
erage of some 13 000 local minima till finding a structure 



(b) 



FIG. 1. The LVPS is illustrated for a Si64 supercell in a per- 
spective (a) and a orthographic (b) view. The yellow (light) 
and red (dark) spheres denote the positions of the atoms be- 
fore and after applying LVPS, respectively. The resulting su- 
percell consists of 66 atoms. The primitive cell found by the 
LVPS is shown in the right bottom corner by a blue paral- 
lelepipede. 



with ground state enthalpy, an improvement in perfor- 
mance by two orders of magnitude compared to the ear- 
lier results. However, since the EDIP potential has only 
first neighbor interaction, the cubic diamond structure 
and its polytypes (for example the hexagonal diamond 
(Lonsdaleite) structure) all give the same ground state 
enthalpy. Therefore, only 12 out of the 45 successful 
runs arrived at the cubic diamond structure. It needs to 
be emphasized that this is not a shortcoming of the MH 
method, but an effect of the short range character of the 
EDIP potential. Increasing the number of distinct and 
accepted minima to be found will lead to an increase in 
the success rate. So, when increasing the search length 
by a factor of 6 the success rate is almost doubled to 80%. 
Another approach to increase the success rate of the 
MH method is by parallelizing the MH runs. In our expe- 
rience, the success rate of a serial MH simulation can de- 
pend heavily on the initial guess of the structure. There- 
fore, it can be advantageous to start multiple parallel 
simulations starting from different initial structures. We 
performed a simple statistical analysis to estimate the 



computational cost necessary if, instead of performing 
multiple serial runs with index i, a parallel run is started 
with m processes of which each requires k MH steps to 
find the ground state. In the parallel version all processes 
would be stopped as soon as one of the processes finds the 
ground state structure. The necessary number of min- 
ima to be visited is given by n™ in = to ■ min{Zi} where 
i = 1, . . . , m. An average value n™ in should be computed 
from runs with different initial structures. Then, the op- 
timal number of processes for the particular problem is 
?7i = 77i resulting in the minimal value of n™ in . Applied 
to the Si64 system the computational cost can be reduced 
to half for too = 6, n^ in rs 6200 

For non-periodic systems the identification of struc- 
tures based on their energies is sufficiently accurate in 
most cases, but the enthalpy degeneracy of the poly- 
types within shortranged potentials requires some addi- 
tional method to distinguish different structures. Since 
the most natural characteristic of a structure is its geome- 
try we use a fingerprint function related to experimental 
diffraction patterns proposed by Valle et al.M- A con- 
tinuous one-dimensional function is defined by summing 
up weighted Gaussian functions centered at all relative 
atomic distances. To reduce numerical errors the fin- 
gerprint function is discretized into to bins, leading to a 
vector of size to uniquely related to the structure. By 
using the angle between two fingerprint vectors a cosine 
distance is defined which then can be used to determine 
the similarities between structures. 
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2000 4000 6000 8000 10000 
Number of visited minima 

FIG. 2. The enthalpies of the local minima visited during 
a MH simulation of a Si64 supercell at zero pressure, start- 
ing from a random input configuration. The enthalpies were 
shifted such that the ground state enthalpy is zero. The in- 
set shows clearly the crossing of the potential barrier before 
arriving at the ground state structure. 



Frequently, when a global optimization run does not 
reach the global minimum for a long time, the simula- 
tion is stuck in an enthalpy funnel with barriers hard to 
overcome. In most cases of our simulation of Si64 at zero 
pressure these funnels are determined by lattice vectors 
which can fit a cubic diamond structure with 62 or 66 sil- 
icon atoms, but not exactly the given 64 atoms. In these 



cases, a diamond structure (or one of its polytypes) fit 
into the cell perfectly with exception of some defective 
subregions where 2 Si atoms are either missing or are 
redundant. For these cases we developed a lattice vec- 
tor prediciton scheme (LVPS) to modify the simulation 
cell by adding or removing atoms such that a perfect 
crystalline structure is recreated. We define a three di- 
mensional scalar function /(r) by summing up Gaussian 
functions with a width a on each atomic position r^ 
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The Gaussian width should be small enough that the 
overlap between neighboring atoms is vanishingly small. 
The autocorrelation function is defined as 
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In Fourier space this convolution is simply a multipli- 
cation of the individual fourier transformed functions of 
/* and /, and according to the Wiener- Chinchin theorem 



H = T{h) = {T{f)Y-Hf) = \Hft 
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where T denotes the Fourier transform and * denotes the 
complex conjugate. A transformation back to real space 
results in the desired autocorrelation function. The au- 
tocorrelation function is then scaled such that the peak 
at the origin is 1, and periodic boundary conditions are 
applied. The three peaks closest to the origin with an am- 
plitude of more than 0.5 spanning a parallelepipede with 
a non- vanishing volume give the vectors of a primitive 
unit cell. Using this cell to span the whole simulation 
cell, the most common basis is identified and the crys- 
talline structure is reproduced (see Fig. [J). Using LVPS 
to identify ground state structures the success rate was 
further increased to almost 95%. The LVPS can be fur- 
ther used to predict the correct system size when the 
number of atoms in the crystal basis is not known in 
advance. 

Fig. [5] shows a typical progress of a MH run for a Si64 
supercell at zero pressure. First, starting from a random 
structure, the enthalpy decreases as new parts of the con- 
figurational enthalpy surface are explored. After visiting 
some 4000 local minima, the system is caught for some 
time in a deep funnel. This funnel corresponds exactly 
to the case discussed above and the perfect crystal could 
be completed applying the LVPS. However, after visiting 
slightly more than 10 200 minima a barrier is crossed and 
the system finally reaches the ground state structure. 

The EDIP potential proved to be somewhat inaccurate 
when predicting the enthalpies of structures at a pressure 
of 16 GPa. Both the well-known /3-tin structure 4 ^— (Si- 
ll) and the Imma structure 48 (Si- XI) were found to be 
metastable in EDIP. The simple hexagonal phased (Si- 
VII) is not even a local minimum on the enthalpy sur- 
face and relaxes to the simple cubic structure. Instead, a 



structure with shifted layers of cubic elements where each 
atom is fourfold coordinated with a bond lenght « 2.4 A 
was predicted as the ground state with an enthalpy of 
—2.923 eV and a volume of 14.378 A 3 per atom, a struc- 
ture not in the fitting database used for EDIP. Finding 
this unexpected crystalline ground state shows the pre- 
dictive power of the MH method for unknown structures. 
However, the novel structure was found to be unstable 
within density functional theory calculations. The sec- 
ond lowest crystalline structure was found to be the bct-5 
structure^, with an enthalpy of —2.865 eV and a volume 
of 15.264 A 3 per atom. 




FIG. 3. The type-I silicon clathrate is shown with the unit 
cell represented by the green box. Two types of cages are 
used to compose the structure, a small pentagonal dodecahe- 
dron (blue) and a larger hexagonal truncated trapezohedron 

(red). ^ 

As a next application we studied the global optimiza- 
tion scheme on a Si46 supercell at zero pressure with vol- 
ume constraints which was realized by introducing a har- 
monic energy term E vo i (second term in equation (|16[1 ) 
in addition to the standard EDIP potential energy: 
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where k is a small scalar value and $7o is the target vol- 
ume. The additional gradient on the cell vectors is given 
by 
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The Si46 was chosen since it can form the unit cell of 
a type-I Si46 clathrate structure^! (see Fig. [3]). Al- 
though clathrates have fourfold coordinated structures 
their overall geometry differ significantly from the cubic 
diamond structure. Composed from polyhedral build- 
ing blocks, a major part of the unit cells remains void, 



resulting in porous crystals. Therefore, the type-I S146 
clathrate within the EDIP potential has a large volume 
per atom ratio of 22.852 A 3 . 

Starting from random input positions 20 MH runs were 
started using the equilibrium volume of the type-I Si46 
clathrate as the target fio and k = 2.0, each process visit- 
ing some 1.5 Mio. minima. Only one run was able to find 
the clathrate ground state unit cell after visiting roughly 
150 000 minima. Nevertheless, finding the clathrate unit 
cell at all is an encouraging result since this particular 
system is a big challenge in global optimization for the 
following reasons. First, there is only one unit cell cor- 
responding to the ground state and it consists of a huge 
basis of 46 atoms with a complex structure. Second, there 
are two main funnels that compete during the search for 
the ground state. On one hand, the system prefers to 
crystallize to the cubic diamond structure, but the void 
areas with the dangling bonds are energetically not favor- 
able. On the other hand there is a tendency of forming 
spherical cage-structures in a porous crystal, but it is sel- 
dom possible to obtain tetrahedral bond angles. These 
are two fundamentally different structures and are sepa- 
rated by a very high potential barrier hard to overcome, 
hence starting from many different random input guesses 
is crucial for a successful run. 
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FIG. 4. (a) A cell of the putative ground state of a BLJ 
crystal consisting of 48 type A and 12 type B atoms. Type A 
atoms are denoted by red (large) spheres, type B atoms are 
blue (small). The detailed structure of the embedded type B 
atoms is shown in subfigure (b). The 7 type A atoms form 
a monocapped trigonal prism^ The symmetry group of this 
isolated molecule is Ci v 



D. Application on binary Lennard-Jones 
mixtures 

As a last application we reinvestigated a much-studied 
BLJ mixture at zero pressure, a system widely accepted 
as benchmark systems. Lately, these mixtures have been 
studied by Middelton et al.M- They found that ordered 
crystalline phases are energetically favored, contrary to 
earlier results indicating a preference for glassy, amor- 
phous structures. The putative ground state structures 
found in this previous work are available on the Cam- 
bridge Cluster Database.— We studied a supercell with 
60 atoms consisting of 80% type A and 20% type B com- 
ponents. In our calculations we use a small modification 
of the well-known Lennard-Jones potential^ (also used in 
Ref. [53|), truncated and shifted using a quadratic func- 
tion such that both the energy and the first derivative 
are continuous at the cutoff distanced The functional 
form of the pair potential is given in equation (|18[) where 
the subindices a and /3 denote the atom types A and B, 
r a fj is the interatomic distance, e Q ^ is the potential well 
depth (not to be confused with the strain tensor) and o~ a p 
corresponds to the distance where the potential vanishes. 
The potential is zero when the radial distance is larger 
than the cutoff r™g . All parameters are identical to the 
ones used in Ref. [53], namely: o~aa — 1-00, (jab = 0.80, 
(Tbb = 0.88, €aa = 1-00, 6ab = 1.50, €bb — 0.50, and 
a cutoff radius of r™| = 2.5o- a p. 
thalpies are given in units of caa 



— 2.5a a p. All energies and en- 




5£H +4 ^ (18) 



We found several thousand structures with enthalpies 
lower than the -7.08 caa per atom of the crystal previ- 
ously found in Ref. 53]. The putative ground state was 
found to have an enthalpy of -7.49 €aa per atom. In fact, 
this value is even lower than any other enthalpy of the 
larger supercells which were investigated by Middelton 
et al.. The periodic cell of the ground state structure is 
shown in Fig. |4] (a) . It consists of two different regions 
similar to a phase separation. The lower section consists 
purely of type A atoms in a slightly distorted hexago- 
nal closed packed structure. The upper half of the cell 
consists of a mixed phase where type B atoms are em- 
bedded in cages consisting of 7 A-type atoms. A detailed 
illustration of this behaviour is shown in Fig. HKb). 



III. CONCLUSIONS 

We have generalized the MH method to periodic sys- 
tems using variable cell shape MD for the escape step. 
To cross low-enthalpy barriers more quickly along the 
MD trajectory the velocities are chosen to point along 
soft mode directions, taking also into account the lattice 
parameters. In several applications we show that low 
enthalpy structures can be predicted and ground state 
configurations can be found efficiently. For silicon sys- 
tems the ground state structure at zero pressure could 
be found orders of magnitude faster compared to the ini- 
tial MH method with fixed cell shape. Using parallel 
runs and the LVPS, the success rate can be significantly 
increased. The LVPS also allows to find unknown crys- 
talline structures without knowing the exact number of 
atoms of the unit cell in advance. For the Si46 super- 
cell the type-I clathrate structure was found using vol- 
ume constraints. Small BLJ mixture cells were used to 
demonstrate the effect of softening on the optimization 
performance. For larger BLJ mixtures a much-studied 
composition was reinvestigated and new putative ground 
state structures were found. Overall, the periodic MH 
method has shown to be a promising approach in struc- 
ture prediction. 
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